
###############################################################################
# class definition
###############################################################################

setClass("PolyMatrix", representation(columns="list"))
setClass("PolyList", representation(elements="list"))

###############################################################################
# create polynomial matrix
# Arguments:
#       dime    vector with the dimensions.
#       elem    an instance of the class PolyList with the elements
#                to fill the matrix.
# Value: an instance of the class PolyMatrix.
###############################################################################

PolyMatrix<-function(dime,elem)
{
    A=new("PolyMatrix")
    n=dime[1]
    m=dime[2]
    for (j in 1:m) # across columns
    {
        col=new("PolyList") # create column
        for (i in 1:n)      # append elements
        {
            ind=(i+(j-1)*n-1)%%length(elem)+1
            col@elements=c(col@elements,list(elem[[ind]]))
        }
        A@columns=c(A@columns,list(col)) # append column

    }
    return(A)

}

###############################################################################
# access elements
# We ensure to be able to get the entries of a polynomial matrix
# with the brackets [].
###############################################################################

setMethod(
    "[",
    "PolyMatrix",
    definition=function(x,i,j) {
    
        if ((length(i)>1 & min(i)<0) | (length(j)>1 & min(j)<0)){
            cat("If there are negative values, then length must be 1.")
            return()
        }
        
        if (length(i)==1 & length(j)==1 & min(min(i),min(j))>0){
            return(x@columns[[j]]@elements[[i]])
        }
        else{
            if (min(i)>0){
                I=i
            }
            else{
                I=setdiff(1:(dim(x)[1]),-i)
            }
    
            if (min(j)>0){
                J=j
            }
            else{
                J=setdiff(1:(dim(x)[2]),-j)
            }
                                         
            n=length(I)
            m=length(J)
            C=zeroPM(n,m)
            for (k in 1:n){
                for (l in 1:m){
                    set(C,k,l)=x[I[k],J[l]]                        
                }  
            }          
            return(C)                       
        }
    }
)

###############################################################################
# setter
# We ensure to be able to set the entries of a polynomial matrix with []<-.
###############################################################################

setGeneric("set<-", function(x, i,j, value) standardGeneric("set<-"))
setReplaceMethod("set","PolyMatrix",
   function(x, i, j, value) {
      x@columns[[j]]@elements[[i]] <- value
      x
   }
)

###############################################################################
# show
# Method to show polynomial matrices.
###############################################################################

setMethod(
    "show",
    "PolyMatrix",
    definition=function(object) 
    {
        n=dim(object)[1]
        m=dim(object)[2]
        for (i in 1:n){   
            for (j in 1:m){
                cat("[",i,",",j,"]=",as.character(object[i,j]),"\t")
            }
            cat("\n")
        }


    }
)

###############################################################################
# get dimension
# Method to get the dimension of a polynomial matrix.
###############################################################################

setMethod(
    "dim",
    "PolyMatrix",
    function (x) {
        ncol=length(x@columns)
        nrow=length(x@columns[[1]]@elements)
    return(c(nrow,ncol))
    }
)

###############################################################################
# sum of PolyLists
# We overload the operator + to be able to sum two lists of the same length.
###############################################################################

setMethod(
    "+",
    signature(e1 = "PolyList", e2 = "PolyList"),
    function (e1, e2) {
        if (length(e1@elements)!=length(e2@elements))
        {
            cat("error: list length mismatch")
            return(NULL)
        }
        else
        {
            sumelem=list()
            for (i in 1:length(e1@elements)){
                sumelem=c(sumelem, list(e1@elements[[i]]+e2@elements[[i]])) 
            }
            sum=new("PolyList")
            sum@elements=sumelem
            return(sum)
        }
    }
)

###############################################################################
# sum of PolyMatrices
# We overload the operator + to be able to sum two polymatrices of the same dim.
###############################################################################

setMethod(
    "+",
    signature(e1 = "PolyMatrix", e2 = "PolyMatrix"),
    function (e1, e2) {
        
        nr=dim(A)[1]
        nc=dim(A)[2]    

        sumcols=list()
        for (j in 1:nc){
            sumcols=c(sumcols, e1@columns[[j]]+e2@columns[[j]]) 
        }
        sum=new("PolyMatrix")
        sum@columns=sumcols
        return(sum)     
    }
)

###############################################################################
# product of PolyMatrices
# We overload the operator * to be able to multiply two polymatrices of 
# the right dimensions.
###############################################################################

setMethod(
    "*",
    signature(e1 = "PolyMatrix", e2 = "PolyMatrix"),
    function (e1, e2) {
        
        n=dim(e1)[1]
        m=dim(e1)[2]     
        p=dim(e2)[1]
        q=dim(e2)[2]

        if (m!=p){
            cat("error: number of columns of e1 different form number of rows of e2")
            return()
        }

        C=zeroPM(n,q)

        for (i in 1:n){
            for (j in 1:q){
            
                sum=polynomial(coef=c(0))
                for (k in 1:m){                    
                    sum=sum+e1[i,k]*e2[k,j]
                }
               set(C,i,j)<-sum
    
            }
        }
        return(C)       
    }
)

###############################################################################
# product of scalar by PolyMatrix
# We overload the operator * to be able to multiply a scalar by a PolyMatrix.
###############################################################################

setMethod(
    "*",
    signature(e1 = "numeric", e2 = "PolyMatrix"),
    function (e1, e2) {
        
        n=dim(e2)[1]
        m=dim(e2)[2]      
        cat(n,m)

        B=PolyMatrix(c(n,m),polynomial(coef=c(0)))
 
        for (i in 1:n)
        {

            for (j in 1:m)
            {

                set(B,i,j)<-e1*e2[i,j]
            }
        }
        return(B)       
    }
    
)

###############################################################################
# Power of a Matrix (only for integer scalars)
# We overload the operator ^ to be able to calculate the power of a PolyMatrix.
###############################################################################

setMethod(
    "^",
    signature(e1 = "PolyMatrix", e2 = "numeric"),
    function (e1, e2) {
            
        B=e1
        
        for (i in 2:as.integer(e2)){
            B=B*e1   
        }
        return(B)       
    }
)

###############################################################################
# zeroPM Zero Matrix
# Arguments:
#       n    number of rows.
#       m    number of columns.
# Value: an instance of the class PolyMatrix with all entries equal to zero.
###############################################################################

zeroPM<-function(n,m)
{
    return(PolyMatrix(c(n,m),polynomial(c(0))))
}

###############################################################################
# unitPM Unit matrix
# Arguments:
#       n    number of rows and columns.
#
# Value: an n x n unit PolyMatrix.
###############################################################################

unitPM<-function(n)
{
    p0=polynomial(c(0))
    p1=polynomial(c(1))
    listpol1=list(p1)
    for (i in 1:n)
    {
        listpol1=c(listpol1,list(p0))
    }
    
    listpol2=list()

    for (i in 1:(n-1))
    {
        listpol2=c(listpol2,listpol1)
    }


    return(PolyMatrix(c(n,n),listpol2))
}
###############################################################################
# transposition
# We overload the method t to be able to calculate the transpose
#     of a PolyMatrix.
###############################################################################

setMethod(
    "t",
    "PolyMatrix",
    definition=function(x) 
    {
        B=x
        for (i in 1:length(x@columns))
            for (j in 1:length(x@columns[[i]]@elements))
            {
                B@columns[[i]]@elements[[j]]=x@columns[[j]]@elements[[i]]
            }
                
        return(B)

    }
)

###############################################################################
# determinant
# We overload the method det to be able to calculate the
#    determinant of a PolyMatrix.
###############################################################################

setMethod(
    "det",
    "PolyMatrix",
    definition=function(x){        
        n=dim(x)[1]
        m=dim(x)[2]
        if (n!=m){
            cat("Determinant can only be computed for square matrices\n")
            return()
        }
        if (n==1){
            return(x[1,1])        
        }
        else{
            pol=polynomial(coef=c(0))
            for (i in 1:n){        
                pol=pol+((-1)^(i+1))*x[i,1]*det(x[-i,-1])   
            }
        }
    return(pol)
    }
)

###############################################################################
# trimpol We purge a polynomial from higher-order coeffcients that are smaller 
#           than a certain tolerance.
# Arguments:
#       pol    polynomial.
#       tol    tolerance.
#
# Value: a polynomial with equal or less degree.
###############################################################################

trimpol<-function(pol,tol){
    c=coef(pol)
    z=(abs(c)<tol)
    if (all(z)){
        deg=0
    }
    else{
        deg=max(which(!z))-1
    }
    c[z]=0
    return(polynomial(coef=c[1:(deg+1)])) 
}

###############################################################################
# trimpol We purge the entries of a PolyMatrix.
# Arguments:
#       A      a PolyMatrix.
#       tol    tolerance.
#
# Value: a PolyMatrix with entries of equal or less degree.
###############################################################################

trimPolyMat<-function(A,tol){
    B=A
    n=dim(A)[1]
    m=dim(A)[2]
    
    for (i in 1:n)
        for (j in 1:m){
            set(B,i,j)<-trimpol(A[i,j],tol)
            }
            
    return(B)
}

###############################################################################
# elemOperA Elementary operation matrix. If multiplied from the left,
#    adds q times row j to row i. From the right and transposed it performs
#    the same operation by columns.
#
# Arguments:
#       n      dimension.
#       i      row to be changed.  
#       j      row to be added to row i times q.
#       q      scalar.   
#
# Value: a PolyMatrix with units in the main diagonal and q in (i,j).
###############################################################################

elemOperA<-function(n,i,j,q){
    A=unitPM(n)
    set(A,i,j)<-q
    return(A)
    
}

###############################################################################
# elemOperS Elementary operation matrix. If multiplied from the left,
#    switches rows i and j. From the right and transposed it performs
#    the same operation by columns.
#
# Arguments:
#       n      dimension.
#       i      row to be changed.  
#       j      row to be added to row i times q.
#
# Value: a unit PolyMatrix with rows and columns i and j reversed.
###############################################################################

elemOperS<-function(n,i,j){
    A=unitPM(n)
    set(A,i,j)<-polynomial(coef=c(1))
    set(A,j,i)<-polynomial(coef=c(1))
    set(A,i,i)<-polynomial(coef=c(0))
    set(A,j,j)<-polynomial(coef=c(0))
    return(A)
    
}

###############################################################################
# abs_r Size of a polynomial, weighted according to a harmonic sequence.
#
# Arguments:
#       p      polynomial.
#
# Value: numeric value of the size.
###############################################################################

abs_r<-function(p){
    n=length(p)
    w=1/(1:n)
    return(sum(w*abs(as.vector(p))))
    #return(max(abs(as.vector(p))))
}

###############################################################################
# smith Estimate Smith form with the approximate algorithm.
#
# Arguments:
#       A      PoyMatrix.
#
# Value: list(U=U,D=D,V=V), where D is the estimated Smith form and U and V
# are the transforming matrices, so that A=U*D*V
###############################################################################

smith= function(A,tol,div) 
{
    rho = 1 + tol 
    tol2 = 1e-10 
    n = dim(A)[1]
    m = dim(A)[2] 
    U = unitPM(m)
    V = unitPM(n) 

    cont <- 0 
    maxoper <- 500 
    for (i in 1:(n-1)) { 
        flag_ii_div = T
        while (flag_ii_div)  { 
            flag_ii_z = T 
            while (flag_ii_z)  { 
                flag_row = T 
                while (flag_row)  { 
                    nzcount = 0 
                    for (j in (i+1):n) {
                        if (abs_r(A[i,j]) > tol) { 
                            nzcount = nzcount + 1 
                            q=A[i,j]/A[i,i]
                            r=A[i,j]%%A[i,i]

                            if (abs_r(r) < tol) { 
                            
                                T1=elemOperA(n,i,j,-q)
                                T1inv=elemOperA(n,i,j,q)
                                A=trimPolyMat(A,tol2)
                          
                                A_ant = A 
                                A = A*T1
                                U = T1inv*U
                                cont <- cont + 1 
                                
                            if (cont > maxoper)  { 
                                stop('smith_polymath:MaxOperExceeded','Max Oper Exceeded') 
                            } 
                        } else { 
                            T1 =elemOperA(n,i,j,-q)
                            T1inv=elemOperA(n,i,j,q)
                            T2=elemOperS(n,i,j)
                            A_ant = A 
                            A = A*T1*T2
                            A=trimPolyMat(A,tol2) 
                            U = t(T2)*T1inv*U 
                            cont = cont + 1 
                            if (cont > maxoper)  { 
                                stop('smith_polymath:MaxOperExceeded','Max Oper Exceeded') 
                            } 
                        } 
                    } 
                } 
                if (nzcount == 0)  { 
                    flag_row = F
                } 
            } 
            flag_col= T
            while (flag_col)  { 
                nzcount = 0 
                for (j in (i+1):n) {
                    if (abs_r(A[j,i]) > tol)  { 
                    
                        nzcount = nzcount + 1 
                        q=A[j,i]/A[i,i]
                        r=A[j,i]%%A[i,i]

                        if (abs_r(r) < tol)  { 
                            T1 =elemOperA(n,j,i,-q)
                            T1inv=elemOperA(n,j,i,q)
                            A_ant = A 
                            A = T1*A
                            A=trimPolyMat(A,tol2) 
                            V = V*T1inv
                            cont = cont + 1 
                            if (cont > maxoper)  { 
                                stop('smith_polymath:MaxOperExceeded','Max Oper Exceeded') 
                            } 
                        } else { 
                            T1 =elemOperA(n,j,i,-q)
                            T1inv=elemOperA(n,j,i,q)
                            T2=elemOperS(n,i,j)
                            A_ant = A 
                            A =T2*T1*A
                            A=trimPolyMat(A,tol2)
                            V = V*T1inv*t(T2) 
                            cont = cont + 1 
                            if (cont > maxoper)  { 
                                stop('smith_polymath:MaxOperExceeded','Max Oper Exceeded') 
                            } 
                        } 
                    } 
                } 
                if (nzcount == 0)  { 
                    flag_col = F 
                } 
            } 
            flag_ii_z = F
            for (j in (i+1):n) {
                if (abs_r(A[i,j]) > tol | abs_r(A[j,i]) > tol)  { 
                    flag_ii_z = T 
                } 
            } 
        } 

        flag_ii_div = F
	  if (div){	 
           for (j in (i+1):n) {
             for (k in (i+1):n) {
                q=A[j,k]/A[i,i]
                r=A[j,k]%%A[i,i]
                if (abs_r(r) > tol)  { 
                    flag_ii_div = T
                    T1=elemOperA(n,i,j,polynomial(coef=c(1)))
                    T1inv=elemOperA(n,i,j,polynomial(coef=c(-11)))
                    
                    A_ant = A 
                    A = T1*A
                    V = V*T1inv
                    cont = cont + 1 
                    if (cont > maxoper)  { 
                        stop('smith_polymath:MaxOperExceeded','Max Oper Exceeded') 
                    } 
                    break 
                } 
             } 
             if (flag_ii_div)  { 
                 break 
             } 
           } 
	  }
    } 
} 
D <- A 
return(list(U=U,D=D,V=V))
} 

###############################################################################
# maproots Projects a vector of roots onto the set
# {exp(2*pi*i*k/s): k=0, ..., s-1}.
#
# Arguments:
#       r      vector of roots
#       s      seasonality period .
#
# Value: vector of projected roots.
###############################################################################

maproots<-function(r,s){

    if (length(r)==0){
        return()
    }
    else
    {
        nr=length(r)
        rv=exp(2*1i*pi*(0:(s-1))/s)
        dif=abs(r%*%t(rep(1,s))-rep(1,nr)%*%t(rv))
        assig=apply(dif,1,"which.min")
        return(rv[assig])
    }
    
}

###############################################################################
# maproots.mat distributees a vector of roots r in the right positions
# according to PolyMatrix A. Details explained in the article.
#
# Arguments:
#       r      vector of roots
#       A      PolyMatrix.
#
# Value: diagonal PolyMatrix with the roots of r.
###############################################################################

maproots.mat<-function(A,r){

	n<-nrow(A)
	m<-length(r)

	B<-unitPM(n)

	rA<-list()
    rB<-list()
    for (k in (1:n)){
        rB[[k]]=vector(mode="numeric",length=0)
        rA[[k]]<-solve(A[k,k])
    }
    
    mind<-rep(0,n)
    minind<-rep(0,n)
       
    
	for (j in (1:m)){
		for (k in (1:n)){

			if (length(rA[[k]])!=0){
				dif<-abs(r[j]-rA[[k]])
				mind[k]<-min(dif)
				minind[k]<-which.min(dif)
			}
			else{
				mind[k]<-Inf
				minind[k]<-NaN
			}
			

		}
		mink<-which.min(mind)
		rjk<-which.min(abs(r[j]-rA[[mink]]))

        rB[[mink]]<-c(rB[[mink]],r[j])
        rA[[mink]]<-rA[[mink]][setdiff(1:length(rA[[mink]]),rjk)]

	}

    for (k in (1:n)){
        set(B,k,k)<-poly.from.roots(rB[[k]])
    }

    return(B)
}

###############################################################################
# ar2PM Transforms the autoregressive coefficient matrix into a PolyMatrix.
#
# Arguments:
#       model  a list with a component named "ar" that contains the mmatrix.
#
# Value: autoregressive PolyMatrix.
###############################################################################

ar2PM<-function(model){

    p=dim(model$ar)[1]
    n=dim(model$ar)[2]
    
    PHI=zeroPM(n,n)
    
    d=diag(n)
    
    for (i in 1:n){
        for (j in 1:n){
            set(PHI,i,j)<-polynomial(c(d[i,j],-model$ar[,i,j]))
        }
    }
        
    return(PHI)
    
}

###############################################################################
# dfilter Applies filter expressed as a polynomial to a time series.
#
# Arguments:
#       x       time series.
#       delta   filter polynomial.
#
# Value: filtered series.
###############################################################################

dfilter<-function(x,delta){

    y=filter(x,delta,method="convolution",sides=1)

}                      

###############################################################################
# multp Gets multiplicities of the unit roots in a polynomial.
#
# Arguments:
#       p       polynomial.
#       s       seasonality period.
#
# Value: vector with the multiplicities of the roots.
###############################################################################

multp<-function(p,s){       
                 
    r=solve(p)
    
    if (length(r)==0){
        return(rep(0,s/2+1))
    }
    
    nr=length(r) 
    rv=exp(2*1i*pi*(0:(s/2))/s) 
    dif=abs(r%*%t(rep(1,s/2+1))-rep(1,nr)%*%t(rv)) 
    return(apply(abs(dif)<1e-7,2,"sum"))

}              

###############################################################################
# poly.from.mult Calculates a polynomial with unit roots where the
#    multiplicities are given by a vector.
#
# Arguments:
#       m       multiplicities vector.
#       s       seasonality period.
#
# Value: polynomial with the multiplicities according to m.     
###############################################################################

poly.from.mult<-function(m,s){

        rv=exp(2*1i*pi*(0:(s-1))/s)
        rv[abs(Im(rv))<1e-5]<-Re(rv[abs(Im(rv))<1e-5])
                                          
        aux=polynomial(coef=c(1))
        z= polynomial(coef=c(0,1))  
        
        for (i in 1:(s/2+1)){
             if (abs(Im(rv[i]))>1e-5){                 
                aux=aux*((1-2*Re(rv[i])*z+z^2)^m[i])
            }
            else{
                aux=aux*((1-rv[i]*z)^m[i])
            }
        }    
        c=coef(aux)           
        c[which(abs(c)<1e-7)]=0
        
        return(polynomial(coef=c))

}

###############################################################################
# autoGVEC performs the full process of identification and estimation 
#       of the GVEC model.
# 
# Arguments:
#   x            A multivariate series whose model we want to identify.
#   var.order    Order of the VAR model fitted. If equal to "AIC" or "BIC",
#                the order is identified by an information criterion.
#   pvarmax      Max. order of the VAR. Used only if var.order is "AIC" or "BIC".
#   vec.order    Order of the GVEC model fitted. If vec.order is "AIC" or "BIC",
#                    the order is identified by an information criterion.
#   pvecmax      Max. order of GVEC. Used only if vec.order os "AIC" or "BIC".
#   logtrans     If TRUE, log transform applied.
#   eps1         Margin or tolerance for the approximate Smith algorithm.
#   eps2         Margin or tolerance to separate unit roots
#                    (usually, equal to eps1).
#   d            Deterministic effects. It may be a list with a subset
#                    among "i" (intercept), "t" (trend)
#                or equal to "AIC" or "BIC" to use information criteria instead.
#   method       if equal to 0, first version of the algorithm, otherwise,
#                    second version is applied.
#
#   Value.
#    A list with the following:
#   deltas       is a list with the distinct polynomials of the Smith form.
#   pvec         finally chosen order of the GVEC.
#   estcoef      estimated coefficients of the GVEC.
#   pvar         finally chosen order of the VAR.
#   D            Smith form.
###############################################################################

autoGVEC<-function(x,var.order='BIC',pvarmax=2*frequency(y),vec.order='AIC',pvecmax=2*frequency(y),logtrans=T,eps1=log(nrow(x))/sqrt(nrow(x)),eps2=log(nrow(x))/sqrt(nrow(x)),d="BIC",method=0){
                              
    s=frequency(x)       
    k=ncol(x)
    n=nrow(x)  
    
    if (logtrans){
        y=log(x)
    }
    else{
        y=x
    }
    
    # Determination of the VAR order 
    
    if (class(var.order)=="numeric"){
        armod=ar(y,aic=F,order.max=var.order,method="ols")
        pvar=armod$order               
        cat("fijo")                                                                          
    }
    else{
    
        if (var.order=='BIC'){
            vec.BIC=rep(0,pvarmax)
            for (j in 1:pvarmax){
                armod=ar(y,aic=F,order.max=j,method="ols")
                vec.BIC[j]=log(det(armod$var.pred))+k^2*j*log(n)/n 
            }
            pvar=which.min(vec.BIC)
   
        }
        
        if (var.order=='HQ'){
            vec.HQ=rep(0,pvarmax)
            for (j in 1:pvarmax){
                armod=ar(y,aic=F,order.max=j,method="ols")
                vec.HQ[j]=log(det(armod$var.pred))+2*k^2*j*log(log(n))/n 
            }
            pvar=which.min(vec.HQ)
        }
               
        if (var.order=='AIC'){
            armod=ar(y,aic=T,order.max=pvarmax,method="ols")
            pvar=armod$order 
        }    
        
    }
    
    
    cat("pvar=",pvar,"\n")
    show(armod)
    
    # Smith Form

    flag=T               
    cont=0
    
    if (method==0){    
        ident=GVECid(y,pvar,eps1,eps2)
    }
    else{
        ident=GVECid2(y,pvar,eps1,eps2)
    }

    if (class(vec.order)=="numeric"){
        fitmodel=GVECd(y,ident$deltas,vec.order,s,d)                                                                                           
    }
    else{
        vec.AIC=rep(0,pvecmax+1)
        vec.BIC=rep(0,pvecmax+1)
        for (j in 0:pvecmax){
            fitmodel=GVECd(y,ident$deltas,j,s,d)
            vec.AIC[j+1]=log(det(fitmodel$rescov))+k^2*j/n 
            vec.BIC[j+1]=log(det(fitmodel$rescov))+k^2*j*log(n)/n           
        }
        
        
        if (is.nan(vec.BIC[pvecmax+1])){
          qwerty<-1
        }
        if (vec.order=='BIC'){
            p=which.min(vec.BIC)-1
            fitmodel=GVECd(y,ident$deltas,p,s,d)
            
        }
        else{
            p=which.min(vec.AIC)-1
            fitmodel=GVECd(y,ident$deltas,p,s,d)
        }
    }

    GVECmodel=list(deltas=ident$deltas,pvec=p,estcoef=fitmodel$estcoef,pvar=pvar,PHI=ident$PHI,D=ident$D,sigma=fitmodel$rescov,determ=fitmodel$determ,list.Delta=fitmodel$list.Delta)
    
    cat("V2.0 (26/02/2013)")
                        
    return(GVECmodel)

}

###############################################################################
# GVECd fits GVEC models with different deterministic terms and the same lag structure.
# 
# Arguments:
#   y            The multivariate series to fit.
#   deltas       The list with the distinct polynomials of the Smith from identification step..
#   p            Order of the GVEC.
#   s            Seasonality period (4 for quarterly data, 12 for monthly, ...)
#   d            Deterministic Deterministic effects. It may be a list with a subset among "i" (intercept), "t" (trend)
#                   or equal to "AIC" or "BIC" to use information criteria instead.
#
#   Value. A list with the following:
#   estcoef     Estimated coefficients.
#   covcoef     Covariance Matrix of the coefficients (not to be used).
#   res         Residuals.
#   rescov      Covariance matrix of the residuals.
#   determ      Deterministic effects.
#   list.Delta  Polynomials in the RHS of the model.
###############################################################################

GVECd<-function(y,deltas,p,s,d=list()){
    
    k=ncol(y)
    n=nrow(y)

    determ=list("i","t")

    if (class(d)=="list"){
        fitmodel=GVECfit(y,deltas,p,s,d) 
    }
    else{
        if (class(d)=="character"){
            if (d=="AIC"){
                combinat=as.matrix(expand.grid(c(0,1),c(0,2)))
                vec.BIC=rep(0,nrow(combinat))
                for (j in 1:nrow(combinat)){
                    fitmodel=GVECfit(y,deltas,p,s,determ[combinat[j,]])
                    vec.BIC[j]=log(det(fitmodel$rescov))+k^2*j*log(n)/n                
                }
                comb=determ[which.min(vec.BIC)]
                fitmodel=GVECfit(y,deltas,p,s,comb)
                return(fitmodel)
            }
            else{
                combinat=as.matrix(expand.grid(c(0,1),c(0,2)))
                vec.AIC=rep(0,nrow(combinat))
                for (j in 1:nrow(combinat)){
                    fitmodel=GVECfit(y,deltas,p,s,determ[combinat[j,]])
                    vec.AIC[j]=log(det(fitmodel$rescov))+k^2*j*log(n)/n                
                }
                comb=determ[which.min(vec.AIC)]
                fitmodel=GVECfit(y,deltas,p,s,comb)
                return(fitmodel)            
            
            }    
        }
        else{
            cat("error: d must be either a list of deterministic terms or an information criterion (BIC or AIC)")
            return()
        }
    }
}

###############################################################################
# GVECid identifies the GVEC model with the first version of the algorithm
# Arguments:
#   x            A multivariate series whose model we want to identify.
#   p            Order of the VAR model fitted. 
#   eps1         Margin or tolerance for the approximate Smith algorithm.
#   eps2         Margin or tolerance to separate unit roots (usually, equal to eps1).

#   Value.
#    A list with the following:
#   PHI          Estimated autoregressive polynomial matrix.
#   deltas       is a list with the distinct polynomials of the Smith form.
#   D            Smith form.
###############################################################################

GVECid<-function(x,p,eps1,eps2){

    T=length(x)
    n=ncol(x)
    s=frequency(x)
    
    modfit=ar(x, aic = F, order.max = p,method="ols")

    PHI=ar2PM(modfit)

    smithPHI=smith(PHI,eps1,T)
       
    deltas=list()
    
    for (i in 1:n){
    
        allroots=solve(smithPHI$D[i,i])
        unitr=which(abs(allroots)<(1+eps2))
        
        newdelta=poly.from.roots(maproots(allroots[unitr],s))
       
        if (i==1)
            deltas=c(deltas,list(newdelta)) 
        else{
            if (!(identical(deltas[[length(deltas)]],newdelta))){
                deltas=c(deltas,list(newdelta))
            }
                              
        }                    
    }        
    
    return(list(PHI=PHI,deltas=deltas,D=smithPHI$D))
    
}

###############################################################################
# GVECid2 identifies the GVEC model with the second version of the algorithm
# Arguments:
#   x            A multivariate series whose model we want to identify.
#   p            Order of the VAR model fitted. 
#   eps1         Margin or tolerance for the approximate Smith algorithm.
#   eps2         Margin or tolerance to separate unit roots (usually, equal to eps1).
#
#   Value.
#    A list with the following:
#   PHI          Estimated autoregressive polynomial matrix.
#   deltas       is a list with the distinct polynomials of the Smith form.
#   D            Smith form.
###############################################################################

GVECid2<-function(x,p,eps1,eps2){

    T=length(x)
    n=ncol(x)
    s=frequency(x)
    
    modfit<-ar(x, aic = F, order.max = p,method="ols")

    PHI<-ar2PM(modfit)

    # Get roots or VAR

    rv<-solve(det(PHI))

    # Separate unit roots of VAR

    urv<-rv[abs(rv)<1+eps2]

    # Obtain pseudo-Smith form (i.e., where the elements do not divide each other)

    smithPHI0<-smith(PHI,eps1,F)

    # Assign unit roots of VAR to elements

    urvmap<-maproots(urv,s) # first, project then on exp(2*pi*k)

    D1<-maproots.mat(smithPHI0$D,urvmap)

    # Obtain final Smith form

    smithPHI<-smith(D1,0.001,T)
       
    deltas<-list()
    
    for (i in 1:n){
            
        newdelta<-smithPHI$D[i,i]
       
        if (i==1)
            deltas<-c(deltas,list(newdelta)) 
        else{
            if (!(identical(deltas[[length(deltas)]],newdelta))){
                deltas<-c(deltas,list(newdelta))
            }
                              
        }                    
    }        
    
    return(list(PHI=PHI,deltas=deltas,D=smithPHI$D))
    
}

###############################################################################
# GVECfit fits GVEC models with fixed deterministic terms and lag structure.
# 
# Arguments:
#   y            The multivariate series to fit.
#   deltas       The list with the distinct polynomials of the Smith from the
#                    identification step.
#   p            Order of the GVEC.
#   s            Seasonality period (4 for quarterly data, 12 for monthly, ...)
#   d            Deterministic Deterministic effects. It may be a list with a
#                subset among "i" (intercept), "t" (trend) or equal to
#                    "AIC" or "BIC" to use information criteria instead.
#
#   Value. A list with the following:
#   estcoef     Estimated coefficients.
#   covcoef     Covariance Matrix of the coefficients (not to be used).
#   res         Residuals.
#   rescov      Covariance matrix of the residuals.
#   determ      Deterministic effects.
#   list.Delta  Polynomials in the RHS of the model.
###############################################################################

GVECfit<-function(y,deltas,p,s,d=list()){

    T=nrow(y)
    n=ncol(y)
    r=length(deltas)       
    
    sroots=exp(2*1i*pi*(0:(s/2))/s) 
    
    Tr=T-p-length(deltas[[r]])+1       
    
    yr=dfilter(y,deltas[[r]])
          
    if (p>0){

        X<-lag(yr,k=-1)
                
        if (p>1){
            for (i in (2:p)){
                ylag<-lag(yr,k=-i)
                X<-cbind(X,ylag)
            }
        }
    }
    else{
        X=numeric()
    }                    
    
    if (r>1){
        rightside<-(2:r)
    }
    else
    {
        rightside<-numeric()
    }
    
    contdiff<-0
    list.Delta<-list()
                    
    for (i in rightside){         

        quot= deltas[[i]]/deltas[[i-1]]    
                       
        mq=multp(quot,s)   

        for (j in 1:(s/2+1)){
                         
            if (mq[j]>0){                       
                mseq=(1:mq[j])
            }
            else{
                mseq=numeric()
            }
        
            for (k in mseq){    
              
                mfact<-rep(0,s/2+1)   
                mfact[j]<-k 
                    
                m_delta_succ<-multp(deltas[[i]],s)    
                                        
                m_delta_ijk<-m_delta_succ-mfact       
                
                delta_ijk<-poly.from.mult(m_delta_ijk,s)*polynomial(coef=c(0,1))   
                contdiff<-contdiff+1             
                list.Delta<-c(list.Delta,list(delta_ijk))
                             
                yd<-dfilter(y,delta_ijk)
                        
                if (length(X)>0){
                    X<-cbind(X,yd)                                           
                }
                else
                {
                    X<-yd
                }    
                    
                if (Im(sroots[j])>1e-5){
                                    
                    delta_ijk<-delta_ijk*polynomial(coef=c(0,1))
                    contdiff<-contdiff+1
                    list.Delta<-c(list.Delta,list(delta_ijk))
                    yd<-dfilter(y,delta_ijk)

                    X<-cbind(X,yd) 
                }
            }               
            
        }    
    }            
                   
    lostobs=length(deltas[[r]])+p
    
    #tfirst=c((lostobs-1)%/%s+1,(lostobs-1)%%s+1)
    tfirst=lostobs+1
    tlast=end(yr)
    
    if (length(X)==0){
        return(list(res=yr,rescov=cov(yr)))
    }
    
    Xmat<-window(X,start=tfirst,end=tlast)
    Y<-window(yr,start=tfirst,end=tlast)
    
    # Column names

    for (i in seq(length.out=p))
        for (j in (1:n))
            colnames(Xmat)[n*(i-1)+j]<-paste('ylag',toString(i),'.',toString(j),sep="")

    for (i in seq(length.out=contdiff))
        for (j in (1:n))
            colnames(Xmat)[p*n+n*(i-1)+j]<-paste('yd',toString(i),'.',toString(j),sep="")
    

    if (length(d)!=0){
        
        if ("i" %in% d){
            intercept=ts(data=rep(1,nrow(Xmat)),frequency=s,start=start(Xmat))
            Xmat=cbind(Xmat,intercept)            
        }
        
        if ("t" %in% d){
            trend=ts(data=seq(1,nrow(Xmat)),frequency=s,start=start(Xmat))
            Xmat=cbind(Xmat,trend)            
        }
    
    }    
    
    
    if (any(is.na(Xmat))){
      qwerty<-1
    }

    estcoef=solve(t(Xmat)%*%Xmat)%*%t(Xmat)%*%Y               
    
    res=Y-Xmat%*%estcoef
    Sigma=t(res)%*%res/(T-lostobs+1)
    
  #  covcoef=solve(t(Xmat)%*%Xmat)%*%t(Xmat)%*%(diag(T-lostobs+1)%x%Sigma)%*%Xmat%*%solve(t(Xmat)%*%Xmat)        
  
    covcoef=list()
  
    return(list(estcoef=estcoef,covcoef=covcoef,res=res,rescov=Sigma,determ=d,list.Delta=list.Delta))
    
  
}

###############################################################################
# roots represents the roots of a model.
# Arguments:
#   model        A model list as is returned by autoGVEC.
###############################################################################

roots<-function(model){

    r=solve(det(model$PHI))
    plot(r)
    a <- seq(0,2*pi, length.out=100)
    lines(x=cos(a),y=sin(a))

}

###############################################################################
# GVECpred calculates forecasts.
# Arguments:
#   y            The multivariate series to fit.
#   s            Seasonality period (4 for quarterly data, 12 for monthly, ...)
#   model        A model list as is returned by autoGVEC.
#   h            Number of forecasts to be calculated.
#
# Value:         Array h x dimension of the time series.
###############################################################################

GVECpred<-function(y,s,model,h){

	n=ncol(y)
	pred<-array(dim=c(h,n))

	for (k in (1:h)){
	
		pred[k,]<-GVECpred.1s(y,s,model)
		
		y=ts(data=rbind(y,pred[k,]),start=start(y),frequency=frequency(y))

	}

	return(pred)
}

###############################################################################
# GVECpred.1s calculates one 1-step forecast.
# Arguments:
#   y            The multivariate series to fit.
#   s            Seasonality period (4 for quarterly data, 12 for monthly, ...)
#   model        A model list as is returned by autoGVEC.
#
# Value:        vector with the components of the 1-step forecast.
###############################################################################

GVECpred.1s<-function(y,s,model){

    deltas=model$deltas
    d=model$determ
    p=model$pvec

    n=ncol(y)
    y=ts(data=rbind(y,rep(0,n)),start=start(y),frequency=frequency(y))
    T=nrow(y)
    
    r=length(deltas)       
    
    sroots=exp(2*1i*pi*(0:(s/2))/s) 
    
    Tr=T-p-length(deltas[[r]])+1       
    
    yr=dfilter(y,deltas[[r]])
            
    if (p>0){

        X=lag(yr,k=-1)
        
        if (p>1){
            for (i in (2:p))
                X=cbind(X,lag(yr,k=-i))
        }
    }
    else{
        X=numeric()
    }                    
    
    if (r>1){
        rightside<-(2:r)
    }
    else
    {
        rightside<-numeric()
    }

    for (i in rightside){         

        quot= deltas[[i]]/deltas[[i-1]]    
                       
        mq=multp(quot,s)   

        for (j in 1:(s/2+1)){
                         
            if (mq[j]>0){                       
                mseq=(1:mq[j])
            }
            else{
                mseq=numeric()
            }
        
            for (k in mseq){    
               
                mfact=rep(0,s/2+1)   
                mfact[j]=k 
                    
                m_delta_succ=multp(deltas[[i]],s)    
                                        
                m_delta_ijk=m_delta_succ-mfact       
                
                delta_ijk=poly.from.mult(m_delta_ijk,s)*polynomial(coef=c(0,1))   
                
                if (length(X)>0){
			
                    X=cbind(X,dfilter(y,delta_ijk)) 
			                                          
                }
                else
                {
  
                    X=dfilter(y,delta_ijk)
                }    
                    
                if (Im(sroots[j])>1e-5){

                    delta_ijk=delta_ijk*polynomial(coef=c(0,1))
                    X=cbind(X,dfilter(y,delta_ijk)) 
                }
            }               
            
        }    
    }            

           
    lostobs=length(deltas[[r]])+p
    
    tfirst=c((lostobs-1)%/%s+1,(lostobs-1)%%s+1)
    tlast=end(yr)
    
    if (length(X)==0){
        return(list(res=yr,rescov=cov(yr)))
    }
    
    Xmat=window(X,start=tfirst,end=tlast)
    Y=window(yr,start=tfirst,end=tlast)

    if (length(d)!=0){
        
        if ("i" %in% d){
            intercept=ts(data=rep(1,nrow(Xmat)),frequency=s,start=start(Xmat))
            Xmat=cbind(Xmat,intercept)            
        }
        
        if ("t" %in% d){
            trend=ts(data=seq(1,nrow(Xmat)),frequency=s,start=start(Xmat))
            Xmat=cbind(Xmat,trend)            
        }
    
    }    


    estcoef=model$estcoef             
    
    res=Y-Xmat%*%estcoef
    
    covcoef=list()

    return(res[Tr,])
    
}
